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Abstract 



We propose a new, unified approach to solving jump-diffusion partial integro- 
— i differential equations (PIDEs) that often appear in mathematical finance. Our method 

consists of the following steps. First, a second-order operator splitting on financial pro- 
cesses (diffusion and jumps) is applied to these PIDEs. To solve the diffusion equation, 
we use standard finite-difference methods, which for multi-dimensional problems could 
also include splitting on various dimensions. For the jump part, we transform the 
jump integral into a pseudo-differential operator. Then for various jump models we 
show how to construct an appropriate first and second order approximation on a grid 
which supersets the grid that we used for the diffusion part. These approximations 



make the scheme to be unconditionally stable in time and preserve positivity of the 
solution which is computed via a matrix exponential. The paper demonstrates that 
the proposed method is computationally efficient, accurate and simple to implement. 

1 Introduction 

Partial integro-differential equations (PIDEs) naturally appear in mathematical finance if 
an underlying stochastic process is assumed to be a combination of diffusion and jumps. 



* Opinions expressed in this paper are those of the author, and do not necessarily reflect the views of 
Numerix LLC. 



A wide class of Levy processes fall into this category. In modern popular models such as 
stochastic volatility or, e.g., hybrid models, jumps could accompany any stochastic factor, 
thus increasing the overall complexity of the problem. For more details about jump- diffusion 
processes, see 



Cont and Tankov (2004), Sato (1999). 



Unsurprisingly, most of these PIDEs cannot be solved in closed form. At the same time, a 
numerical counterpart must be efficient. This is especially important if such a jump- diffusion 
model is used not only for pricing (given the values of the model parameters), but for their 
calibration as well. While the solution of the diffusion part (PDE) has been numerously 
discussed in the literature and various methods were proposed (see, e.g., Andersen and 



Andreasen (2000), Brennan and Schwartz (1978), Cont and Voltchkova (2003), Duffy (2006) 



Hout and Foulon (2010), Tavella and Randall (2000)), little can be found for the jump part, 
which according to the Levy-Khinchine formula is represented by a non-local integral. 

Various methods were suggested to address the construction of an efficient algorithm 
for solving these type of PIDEs (see, e.g., a discussion of these methods and problems 



related to their implementation in Carr and Mayo (2007), Itkin and Carr (2012), Strauss 
(2006[)). In particular, they include a discretization of the PIDE that is implicit in the 



differential terms and explicit in the integral term (Cont and Voltchkova 



iterations for computing the integral equation (d'Halluin et al. (2005, 2004)) and a second 



(2003)), Picard 



order accurate, unconditionally-stable operator splitting (ADI) method that does not require 



an iterative solution of an algebraic equation at each time step (Andersen and Andreasen 



( 2000[)). V arious forms of operator splitting technique were also used for this purpose (Itkin 
and Carr (2012)). In this paper, we will review operator splitting on financial processes in 



more detail. 

Assuming that an efficient discretization of the PIDE in time was properly chosen, the 
remaining problem is a fast computation of the jump integral, as it was observed to be 
relatively expansive. We mention three different approaches to numerical computation of 
this integral jj 

The first approach assumes a direct approximation of the integral on an appropriate grid 
and then applies some standard quadrature method, such as Simpson's rule or Gaussian 
quadrature. This approach is computationally expensive for two reasons. First, usually the 
"jump" grid is not the same as the "diffusion" grid. Therefore, after the integral is computed, 
its values at the jump grid should be re-interpolated to the diffusion grid. Second, the integral 
is defined on an infinite domain, so either the domain has to be truncated or a non-uniform 
grid has to be used. Moreover, the complexity becomes greater if an implicit discretization 
of the integral is used, because it requires the solution of a dense system of linear equations 
of a large size. Therefore, most often an explicit discretization is utilized, which brings some 
constraints on the time steps to guarantee stability of the scheme. 

However, an exponential change of variables reduces the expense of evaluating the integral 
at all points. This change converts the integral term into a correlation integral, which can be 
evaluated at all the grid points simultaneously using a Fast Fourier Transform (FFT). This 
approach has been suggested by many authors (Andersen and Andreasen ( 2000[ ), [Tavella and 



1 For some models it can be computed analytically, so in what follows we do not take these models into 



account. 



Randall (2000), Wilmott (1998)). Still, this could be expensive because a large number of 
FFT nodes may need to be used for better accuracy. Another issue is that using FFT to com- 
pute a product of matrix A and vector requires A to be circulant, while the matrix obtained 
after discretization of the jump integral is not that type. Therefore, a direct (naive) usage 
of FFT for this purpose produces undesirable so-called "wrap-around" errors. A common 
technique to eliminate these errors is to embed A, which is actually a Toeplitz matrix, into 
a circulant matrix. This, in turn, requires doubling the initial vector of unknowns, which 



makes the algorithm slower (Itkin and Carr (2012)). Also as noticed in the latter paper for 
GTSP/CGMY/KoBoL models, the FFT algorithm loses its accuracy when the parameter a 
of the CGMY model tends to zeroJ3 

The second approach to computing the jump integral utilizes an alternative representa- 
tion of this integral in the form of a pseudo-differential operator, which puts the entire PIDE 



in the form of a fractional PDE. This problem was considered in Cartea and del Castillo- 



Negrete (2007) and Itkin and Carr (2012). A recent survey of the existing literature on this 



subject and techniques for computation of the jump integral using the Grunwald-Letnikov 



approximation (which is of the first order in space) is given in Andersen and Lipton 



As it is known from Abu-Saman and Assaf (2007), Meerschaert and Tadjeran (2004 



2012 



2006 



Sousa (2008), Tadjeran et al. (2006), a standard Grunwald-Letnikov approximation leads to 



unconditionally unstable schemes. To improve this, a shifted Grunwald-Letnikov approxi- 
mation was proposed, which allows construction of an unconditionally stable scheme of the 
first order in spacefj However, solving pricing equations to second order in the space variable 
is almost an industry standard, and therefore this method requires further investigation to 
address this demand. 

The third methocn exploits a nice idea first proposed in Carr and Mayo (2007). Carr and 
Mayo found that for some Levy models, the solution of the integral evolutionary equationjj 
is equivalent to the solution of a particular PDE. The problem is then to find a proper 
space-differential operator (kernel) to construct such a PDE. Carr and Mayo demonstrated 
the advantage of this approach for the Merton and Kou models, and showed which parabolic 
equations provide the necessary solution. Later in Itkin and Carr (2012), this idea was 



further generalized to the class of pseudo-parabolic equations as applied to a class of Levy 
processes, known as GTSP/CGMY/KoBoL models. These pseudo-parabolic equations could 
be formally analytically solved via a matrix exponential. Itkin and Carr then discuss a 
numerical method for efficiently computing this matrix exponential. When the parameter 
a of the GTSP/CGMY/KoBoL model is an integer, this method uses a finite-difference 
scheme similar to those used for solving parabolic PDEs, and the matrix of this FD scheme 
is banded. Therefore, in this case, the computation of the jump integral: 



2 The CGMY model in this limit becomes the VG model (Madan and Seneta (1990)). 

3 A second-order approximation could in principle be constructed as well, but this would result in a 
massive calculation for the coefficients. Therefore, this appr oach was not further e laborated on. 

4 We also have to mention one more method proposed by Lipton and Sepp (2009J as applied to the Stern 
model. Though it is not evident how to generalize this method for other models, it provides a very efficient 
computational algorithm for the Stern model. 

5 This equation naturally arises at some step of the splitting procedure, if splitting is organized by sepa- 
rating diffusion from jumps. 



Is provided on the same grid as was constructed for the diffusion (parabolic) PDE. 
Outside of this domain (if ever needed, e.g. for European options), the PIDE grid is 
further extended to an infinite domainHbut no interpolation is required afterwards. 



• Has linear (O(N)) complexity in the number of the grid nodes N, since the results 
(e.g., option prices) are given by solving a linear system of equations with a banded 
matrix. In the case of a real parameter a, Itkin and Carr suggested computing the 
prices using the above algorithm at three values of an integer 5 closest to the given 
real a, and then interpolating using any interpolation of the second order. 

In this paper we use a different flavor of this idea. First, we use an operator-splitting 
method on the financial processes, thus separating the computation of the diffusion part 



from the integral part. Then, similar to Itkin and Carr (2012), we represent the jump 



integral in the form of a pseudo-differential operator. Next we formally solve the obtained 
evolutionary partial pseudo-differential equations via a matrix exponential. We then show 
that the matrix exponential can be efficiently computed for many popular Levy models, and 
that the efficiency of this method is not worse than that of the FFT. The proposed method 
is almost universal, i.e., allows computation of PIDEs for various jump- diffusion models in 
a unified form. We also show that this method is relatively simple to implement. 

The rest of the paper is organized as follows. In section[2]we briefly discuss a general form 
of a backward PIDE for the class of Levy models. In Section [3j we introduce a splitting tech- 
nique for nonlinear operators. In section |1J we present our general approach to the solution 
of the PIDE using a splitting and matrix exponential approach. An explicit construction 
of various FD schemes of the first and second order is presented in the next section. There 
we consider the following jump models: Merton, Kou and GTSP (also known as CGMY 
or KoBoL). The results presented in the last two sections are new, and to the best of our 
knowledge have not been discussed in the literature. Our technique utilizes some results from 
matrix analysis related to definitions of M-matrices, Metzler matrices and eventually expo- 
nentially positive matrices. In section |6j we present numerical examples that demonstrate 
the efficiency and accuracy of the proposed method. The final section concludes. 

2 Levy Models and Backward PIDE 

To avoid uncertainty, let us consider the problem of pricing equity options written on a single 
stock. As we will see, this specification does not cause us to lose any generality, but it makes 
the description more practical. We assume an underlying asset (stock) price St be driven by 
an exponential of a Levy process 

S t = S exp(L t ), 0<t<T, (1) 

where t is time, T is option expiration, So = S t \t=o, L t is the Levy process L = (L t )o< t <T 
with a nonzero Brownian (diffusion) part. Under the pricing measure, L t is given by 

Lt = yt + aWt + Y t , 7 ,aGM, a > 0, (2) 



3 In other words the PIDE grid is a superset of the corresponding PDE grid. 



with Levy triplet (7, a, is), where W t is a standard Brownian motion on < t < T and Y t 
is a pure jump process. 

We consider this process under the pricing measure, and therefore e~^ r ~ q ^St is a martin- 
gale, where r is the interest rate and q is a continuous dividend. This allows us to express 7 



as (Eberlein (2009)) 



2 /• 
1 = r-q- — - / (e x - 1 - xl^^i) i/(da 

where is(dx) is a Levy measure which satisfies 



/ e x is(dx) < 00. 

J\x\>l 



>\x\> 

We leave v{dx) unspecified at this time, because we are open to consider all types of 
jumps including those with finite and infinite variation, and finite and infinite activity. [_] 

To price options written on the underlying process St, we want to derive a PIDE that 
describes time evolution of the European option prices C(x,t), x = \og(S t /So). Using a 
standard martingale approach, or by creating a self-financing portfolio, one can derive the 



corresponding PIDE (Cont and Tankov (2004)) 



„. , dC(x,t) 
rC{x,t)= m +1 * 



-a 



+ 



0C(S,t) 1 2 d 2 C(S,t) 
dx 2 dx 2 



C{x + y,t) - C(x,t) - (e y - 1) 



dC(x,t) 
dx 



for all (x,t) G M. x (0,T), subject to the terminal condition 

C(x,T) = h(x), 



is(dy) (3) 



(4) 



where h(x) is the option payoff, and some boundary conditions which depend on the type 
of the option. The solutions of this PIDE usually belong to the class of viscosity solutions 



QCont and Tankov| p004[ )). 

We now rewrite the integral term using the following idea. It is well known from quantum 
mechanics that a translation (shift) operator in L2 space could be represented as 



T a = exp [a 



so 



d_ 

dx t 

Taf(x) = f(x + a). 



7 We recall that a standard Brownian motion already has paths of infinite variation. Therefore, the Levy 
process in Eq.J2| has infinite variation since it contains a continuous martingale component. However, here 
we refer to the infinite variation that comes from the jumps. 



Therefore, the integral in Eq. pi) can be formally rewritten as 



[C(x + y,t)-C(x,t)-(e»-l)^^ 



J = 







exp|^]-l-(^-U^ 



u(dy) = JC(x,t), 

d 



(5) 



u(dy). 



In the definition of operator J (which is actually an infinitesimal generator of the jump 
process), the integral can be formally computed under some mild assumptions about exis- 
tence and convergence if one treats the term d/dx as a constant. Therefore, operator J can 
be considered as some generalized function of the differential operator d x . We can also treat 
J as a pseudo-differential operator. 

With allowance for this representation, the whole PIDE in the Eq.pl) can be rewritten 
in operator form as 

d T C(x,r) = [V + J]C(x,r), (6) 

where t = T — t and T> represents a differential (parabolic) operator 



V 



r a 

2 



d_ 
dx 



1 2 &_ 
2 a dx 2 ' 



(7) 



where the operator T> is an infinitesimal generator of diffusion. 

Notice that for jumps with finite variation and finite activity, the last two terms in 
the definition of the jump integral J in Eq.(|3]) could be integrated out and added to the 
definition of T>. In the case of jumps with finite variation and infinite activity, the last term 
could be integrated out. However, here we will leave these terms under the integral for 
two reasons: i) this transformation (moving some terms under the integral to the diffusion 
operator) does not affect our method of computation of the integral, and ii) adding these 
terms to the operator T> negatively influences the stability of the finite-difference scheme 
used to solve the parabolic equation VC(x,t) = 0. This equation naturally appears as a 
part of our splitting method, which is discussed in the next section. 



3 Operator Splitting Technique 

To solve Eq. (|6~j) we use splitting. This technique is also known as the method of fractional 



steps (see Dyakonov (1964), Samarski (1964), Yanenko (1971)) and sometimes is cited in 



financial literature as Russian splitting or locally one-dimensionally schemes (LOD) (Duffy 
02006) ). 

The method of fractional steps reduces the solution of the original fc-dimensional unsteady 
problem to the solution of k one-dimensional equations per time step. For example, consider 
a two-dimensional diffusion equation with a solution obtained by using some finite-difference 
method. At every time step, a standard discretization on space variables is applied, such 
that the FD grid contains Aq nodes in the first dimension and A*2 nodes in the second 
dimension. Then the problem is solving a system of Aq x A^ linear equations, and the 



6 



matrix of this system is block-diagonal. In contrast, utilization of splitting results in, e.g., 
Ni systems of N2 linear equations, where the matrix of each system is banded (tridiagonal). 
The latter approach is easy to implement and, more importantly, provides significantly better 
performance. 



The previous procedure uses operator splitting in different dimensions. Marchuk (1975) 



and then Strang (1968) extended this idea for complex physical processes (for instance, 
diffusion in the chemically reacting gas, or the advection-diffusion problem). In addition 
to (or instead of) splitting on spatial coordinates, they also proposed splitting the equation 
into physical processes that differ in nature, for instance, convection and diffusion. This idea 
becomes especially efficient if the characteristic times of evolution (relaxation time) of such 
processes are significantly different. 

For a general approach to splitting techniques using Lie algebras, we refer the reader 
to Lanser and Verwer (1999). Decomposing the total (compound) operator C for problems 



of interest seems natural if, say, C can be represented as a sum of k noncommuting linear 
operators ^ i=1 A- in this case the operator equation Cf{t) = 5^i=iA/(*) = 0, with 
f(t) being the unknown dependent variable, can be formally integrated via an operator 
exponential, i.e., 



/(*) = c"/(0) 



,£ ti At 



/(o). 



Due to the noncommuting property, the latter expression can be factorized into a product 
of operators 

f(t) = e c K..e^f(0). 

This equation can then be solved in N steps sequentially by the following procedure: 

/ (1) = e £l /(0), 
/(2 ) = e r 2/ « 



CO 






C k f (fc-l) 



(k) 



This algorithm is exact (no bias) if all the operators commute. If, however, they do not 
commute, the above algorithm provides only a first-order approximation in time (i.e., 0(t)) 
to the exact solution. 

To get the second-order splitting for noncommuting operators, Strang proposed a new 



scheme, which in the simplest case (k = 2) is (Strang (1968)) 

f(t) = e Lt f(0) = e^ Ll+L ^f{0) « e Ll te L2 V^/(0) + 0(f). 



(8) 



For parabolic equations with constant coefficients, this composite algorithm is second- 
order accurate in t provided the numerical procedure that solves a corresponding equation 
at each splitting step is at least second-order accurate. 



For nonlinear operators, the situation is more delicate. As shown in Koch and Thalham- 



mer 



(2011). the theoretical analysis of the nonlinear initial value problem 

u\t) = F(u(t)), 0<t<T 

for a Banach-space- valued function u : [0,T] — > X given an initial condition u(0) could be 
done using calculus of Lie derivatives. A formal linear representation of the exact solution is 



u{t) =S F {t,u{0)) 



JD f 



«(0), 



< t < T, 



where the evolution operator and Lie derivatives are given by 



e tD ^v 



S F (t, v), e WF Gv = G(£ F (t, v)), < t < T, 



D F v = F(v), D F Gv = G'(v)F(v) 



for an unbounded nonlinear operator G : D(G) Cl->I. Using this formalism, Koch and 



Thalhammer (2011) showed that Strang's second-order splitting method remains unchanged 



in the case of nonlinear operators. 

Using this result for Eq.d6| gives rise to the following numerical scheme: 



C«( 



e 2 



■v 



C(x,t), 



(9) 



C(x,T + Ar)=e^C (2) (x,r). 



Thus, instead of an unsteady PIDE, we obtain one PIDE with no drift and diffusion (the 
second equation in Eq.p])) and two unsteady PDEs (the first and third ones in Eq.p])). 

In what follows, we consider how to efficiently solve the second equation, while assuming 
that the solution of the first and the third equations can be obtained using any finite- 
difference method that is sufficiently efficient. To this end, in various examples given in the 
next sections we will explicitly mention what particular method was used for this purpose. 

In this paper, we do not discuss the uniqueness and existence of the solution for the 
PIDE; to do so would move us to the definition of a viscosity solution for this class of 



integro-differential equations. For more details, see Cont and Tankov (2004) and Arisawa 



(2005) 



Lastly, let us mention that J = (j)(—id x ), where <p{u) is the characteristic exponent of 
the jump process. This directly follows from the Levy-Khinchine theorem. 



4 Solution of a Pure Jump Equation 

We begin with the following observation. By definition of the jump generator J , under 
some mild constraints on its existence, J could be viewed as a function of the operator d x . 
Therefore, solving the integral (second) equation in Eq.(|9]) requires a few steps. 

First, an appropriate discrete grid G(x) has to be constructed in the truncated (originally 
infinite) space domain. This grid could be nonuniform. An important point is that in the 



space domain where the parabolic equations of Eq. ^ are defined, this grid should coincide 
with the finite-difference grid constructed for the solution of these parabolic equationsjj This 
is to avoid interpolation of the solution that is obtained on the jump grid (the second step 
of the splitting algorithm) to the diffusion grid that is constructed to obtain solutions at the 
first and third splitting steps. 

To make this transparent, let the parabolic equation be solved at the space domain 
[xo,Xfc], xq > — oo, Xk < oo using a nonuniform grid with k + 1 nodes (xo,Xi,...,Xk) and 
space steps hi = X\ — Xo, ■•-, hk = Xk — Xk-i- The particular choice of xq and Xk is de- 
termined by the problem under consideration. We certainly want \xq\ and \xk\ not to be 
too large. The integration limits of J in Eq.([5]) are, however, plus and minus infinity. 
Truncation of these limits usually is done to fit memory and performance requirements. 
On the other hand, we want a fine grid close to the option strike for better accuracy. 
Therefore, a reasonable way to construct a jump grid is as follows. For x$ < x < Xk, 
the jump grid coincides with the grid used for solution of the parabolic PDEs. Outside of 
this domain, the grid is expanded by adding nonuniform steps; i.e., the entire jump grid is 
x_K, %i-k, ■■■X-i,Xq, Xi, ..., Xk, Xk+i, •••, Xk+M- Here K > 0, M > are some integer numbers 
that are chosen based on our preferences. Since contribution to J from very large values of 
x is negligible, the outer gridpoints X-k,%i-k, ■■■X-i and Xk+i, ...,Xk+M can be made highly 
nonuniform. One possible algorithm could be to have the steps of these grids be a geometric 
progression. This allows one to cover the truncated infinite interval with a reasonably small 
number of nodes. 

Second, the discretization of d x should be chosen on G(x). We want this discretization 
to: 

1. Provide the necessary order of approximation of the whole operator J in space. 

2. Provide unconditional stability of the solution of the second equation in Eq.([9]). 

3. Provide positivity of the solution. 

Let A x denote a discrete analog of d x obtained by discretization of d x on the grid G(x). 
Accordingly, let us define the matrix J(A X ) to be the discrete analog of the operator J on 
the grid G(x). The following proposition translates the above requirements to the conditions 
on J(A X ). 

Proposition 4.1 The finite- difference scheme 

C(i,r + Ar) = e ArJ ^C(x,r) (10) 

is unconditionally stable in time t and preserves positivity of the vector C(x, r) if there exists 
an M-matrix B such that J(A X ) = —B. 



Proof By definition of an M-matrix (see Berman and Plemmons (1994)), the class of M 



matrices contains those matrices whose off-diagonal entries are less than or equal to zero, 



^So the PIDE grid is a superset of the PDE grid. 



while all diagonal elements are positive. All eigenvalues of an M-matrix have a positive 
real part. Therefore, if B is an M-matrix, all eigenvalues of J(A X ) have a negative real 



part. Therefore, ||e "^ *)|| < 1 (in the spectral norm), and thus the scheme Eq.(lO) is 
unconditionally stable. 



Now since B is an M-matrix, J is a Metzler matrix (Berman and Plemmons (|1994)) 



An exponential function of the Metzler matrix is a positive matrix. Therefore, if C(x, r) is 



positive, the scheme Eq.(lO) preserves the positivity of C(x, r + At) 



This proposition gives us a recipe for the construction of the appropriate discretization 
of the operator J . In the next section, we will give some explicit examples of this approach. 

Once the discretization is performed, all we need is to compute a matrix exponential 
e Arj(A x )^ anc j ^hen a product of this exponential with C(x, r). The following facts make this 
method competitive with those briefly described in the introduction. We need to take into 
account that: 

1. The matrix J(A X ) can be precomputed once the finite-difference grid G(x) has been 
built. 

2. If a constant time step is used for computations, the matrix A = e ArJ ( AxS) can also be 
precomputed. 

3. If the above two statements are true, the second splitting step results in computing a 
product of a matrix with time-independent entries and a vector. The complexity of 
this operation is 0(N 2 ), assuming the matrix A is N x N, and the vector is N x 1. 
However, N in this case is relatively smalljj One can compare this with the FFT 



algorithm proposed in Andersen and Andreasen (2000) to compute the correlation 



integral. This translates into computation of two matrix-by-vector products. This 
algorithm is 2c x 0(Nlog 2 N), where c is some coefficient. However, N is relatively 
high in this case. Typical values are N = 4096. On top of that, as was discussed 
in the introduction to avoid "wrap-around" effects this number must be doubled, i.e., 
iV = 8192. Also a post-solution interpolation is required]^] Finally, for some models 
(CGMY, VG), the computation of the integral in a neighborhood of x = requires 



special treatment (Cont and Voltchkova (2003)). Overall, the total complexity of FFT 



with allowance for all these remarks becomes almost the same, or even worse than that 
of the proposed matrix exponential algorithm. 

The above consideration is sufficiently general in the sense that it covers any particular 
jump model where jumps are modeled as an exponential Levy process. Clearly, for some 
models computation of the jump integral can be readily simplified. For instance, Merton's 
model, which we discuss in the next section, allows another approach with a better perfor- 
mance. Below we discuss this approach in more detail as it seems to be general enough and 
applicable to some other models as well. 



9 Typical values are N = 100 - 200. 

10 In more advanced approaches, this step could be eliminated; see Parrot (2009 1 . 
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5 Examples for Some Popular Models 



In this section, we review some popular jump models known in the financial literature. Given 
a model, our goal is to construct a finite-difference scheme, first for A x , and then for J(A X ), 



that satisfies the conditions of Proposition 4.1 



5.1 Merton Model 



Merton (1976) considered jumps that are normally distributed with the Levy density 

(x- hj) 2 ' 



v(dx) = A 



1 



27TCTJ 



cxp 



2a] 



dx, 



(11) 



where A, fij and <jj are parameters of the model. Considering the pure jump part of the 
Merton model, one can see that it exhibits finite activity, i.e., a finite number of jumps within 
any finite time interval. Plugging Eq.(ll ) into the definition of the operator J in Eq.(J5J) and 
fulfilling a formal integration gives 



J = X (g/^ds + 5 



-hajdl 



kgL 



K 



a j 



(12) 



where d x 
solved is 



d/dx, d\ = d 2 /dx 2 . The corresponding evolutionary pure jump equation to be 



C(2), 



X,T 



AC^{ 



X,T) 



A = exp 



AAr ( e ^+\° 2 A - K d v -l 



(13) 



A matrix exponential method for this model with the exponential operator as in Eq. ( 13 ) 



has already been considered in Tangman et al. (2011) using a different derivation (from Carr 
and Mayo ( 2007[ )). They also discuss in more detail various modern methods for computing 
the matrix exponentials. 

Recall that the diffusion equations in Eq.(J9J) have to be solved up to some order of 
approximation in time r. Suppose for this purpose we want to use a finite-difference scheme 
that provides a second-order approximation, 0((Ar) 2 ). However, Eq.(13) gives an exact 



solution of the corresponding pure jump equation (the second step in Strang's splitting 
scheme). Since Strang's scheme guarantees only second-order accuracy (0((Ar) 2 )) to the 
exact solution of the full PIDE, the second step could be computed to the same order of 
accuracy. 

Suppose first that we require only an O(Ar) approximation. To this end we can use the 
(1,0) Pade approximation of e ArJ , 

e ArJ w 1 + At J + O(Ar), 

which actually is a well-known Euler explicit scheme. Now the product 

JC {1 \x, t) = -\{Kd x + l)C {l) {x, t) + Ae^ + M d "C (1) (x, r ) 



ii 



It is actually a double exponential operator. 
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can be efficiently computed if one observes that vector 

z (x,t) =e" J * B+ M<£c( 1 )(ar,T) 

is a solution of 

dz{x,s) ( 1 2 2 \ 

g s = [Hjdx + 2 a J d x) z\p, 8) 

for < s < 1 and z(x, 0) = C^'(x,t). This PDE is actually a heat equation and can be 
solved with complexity O(N). In other words, we just rederived the result and method first 
obtained in Carr and Mayo ( |2007 ). 

If one needs a second-order approximation in time, one can use a (1,1) Pade approxima- 
tion of e Ar ^ , 



-AtJ 



so that C(x, t) at this step is a solution of 

(l - \^rj\ C^ 2 \x,r) =(l+ \Atj\ CW(x, 



T . 



Since after discretization of J is performed, the discrete matrix J (Ax) in the right-hand 
side of this equation is dense, a straightforward attempt to solve this matrix equation is not 



efficient. However, it could be solved by using Picard iteration. To start the iteration, one 
needs to compute a product J(Ax)C^(x,r). At every iteration i, one computes a product 



J(Ax)C^_ l (x, t), where C\ (x,r) is the ith approximation of C^- 2 \x,r) obtained after i 
iterative steps. Both these products can be computed by solving a kind of a heat equation 
(as was mentioned above) with complexity O(N). 



5.2 Kou Model 



The Kou model, proposed in Kou and Wang (2004), is a double exponential jump model. 
Its Levy density is 



v{dx) = A (p9 1 e- e ' x l x >o + (1 - p)e 2 e e2X l x<0 ) dx, 



(14) 



where 9i > 1, 9 2 > 0, 1 > p > 0; the first condition was imposed to ensure that the stock 
price S(t) has finite expectation. Using this density in the definition of the operator J in 
Eq.([5]) and carrying out the integration (recalling that we treat d/dx as a constant) gives 



J = \ [-\ + ^ a + p6 l (6 l 
p 



p)9 2 (a + 9 2 )- 1 ] 



(15) 



a = d r 



/'o 



-h zr-r^-, -9 2 < Re{a) < 6 X . 



Ox-l l + 6 : 



See, for instance, d'Halluin et al. (2005), where conditions were derived on approximation to J for 



convergence of these iterations. 
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The inequality —9 2 < Re(a) < 9% is an existence condition for the integral defining J and 
should be treated as follows: the discretization of the operator a = d x should be such that 
all eigenvalues of matrix A, a discrete analog of a, obey this condition. 

Below we assume that A > 0. The construction of the finite-difference scheme in the case 
A < can be done in a similar way. 

Let us first consider discretizing a using first-order approximation in hi, i — 1, . . . ,N, 
where {hi} are the step sizes of the nonuniform space grid. Observe that if 9\ I — A is an M- 
matrix^ then its inverse is a positive matrix. Therefore, a one-sided forward discretization 
of a, which we denote as A F : dC/dx = [C(x + h,t) — C(x,t)]/h, should be chosen for 
this term. In contrast, for the term 9 2 I + A to be an M-matrix, a one-sided backward 
discretization of a, denoted as A B , on the same grid needs to be applied. As now both 
Mi = p9i(9il — A F ) _1 and M 2 = (1 — p)9 2 {9 2 I + A B )~ l are positive matrices, we need the 
matrix Mq = (iqA — I (where A is some discretization of a) to be a Metzler matrix in order 
to make the entire sum J(A X ) = X(Mq + Mi + M 2 ) a Metzler matrix. As /io > 0, this means 
that M = fioA F — I. Thus, by construction, we have proved the following proposition: 



Proposition 5.1 Suppose that discretization of Eq.(15) is done using first-order approxi 



mation of a in hi, i = 1, . . . , N. Also suppose that the following scheme is performed: the 



operator a in the second and third terms in the right-hand side of Eq.(15) are approximated 
using one-sided forward differences, and the operator a in the fourth term is approximated 
using a one-sided backward difference. Then the matrix J(A X ) = A(M + Mi+M 2 ), a discrete 
counterpart of the operator J , is a Metzler matrix. 



According to Proposition |4.1[ all that remains to prove is that the diagonal elements of 
J(A X ) are negative, because then J(A X ) is an M-matrix with a negative sign. In doing so, 
first observe that both A F and A B are triangular with nonzero elements on the diagonal. 
For simplicity, let us consider a uniform grid in x with step size h. Also let us represent Mi 
in the form 

M 1 =p9 l h(9 1 hI-U F )-\ 

where the matrix U F is nonzero only on the main diagonal (where all elements are -1) and 
the first upper diagonal (where all elements are 1). Accordingly, we represent M 2 as 

M 2 = (l-p)9 2 h(9 2 hI + U B )-\ 

with U B = (U F ) T . Obviously all eigenvalues of U are -1, and those of U B are 1. Therefore, 
the diagonal elements of Mi are p9\hj{9\h + 1), and those of M 2 are (1 — p)9 2 h{9 2 h + 1). 
Thus, the diagonal elements of J(A X ) are 



J{A x )i } i — A 



Mn „ h , h 

-1 - ^+p9i— + (1 -p)9 2 . 



h r h9i + l y r ' h9 2 + l 



N]. (16) 



Now the last two terms can be made as small as necessary by choosing a sufficiently small 
value of /ipjso that J(A x )i^ < for alH = 1, . . . , N. Thus, by Proposition 4.1 our scheme is 



13 I here denotes a unit matrix. 

14 The upper limit for H can be found by equating the right-hand side of Eq.(16) to zero and solving the 
resulting equation for h. 
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unconditionally stable, preserves positivity of the solution and is of first-order approximation 
in h. 

The same result can be readily obtained in the case of a nonuniform grid. The discretiza- 
tion of a is constructed in a similar way; instead of a constant grid step size h, variable step 
sizes hi should be used. 

Again, after an appropriate approximation of d x is chosen, J(A X ) can be precomputed 
just once. 

The following proposition extends this method to construct a second order approximation 
scheme for the operator J . 

Proposition 5.2 Suppose that the following discretization scheme is performed: 



J(A X ) + XaA 2 = X{ - I + /j A c + aA 2 + 



pOi 



Oil 

1 



+ 



[i- P )e 2 



h 



a f e 2 i + A B 

pOi(A F ) 2 (1_ 

[dj-AF) 2 



p)0 2 (A 



B\2 



ii + A B y 



(17) 



where a > is a constant, A c is a centered finite- difference approximation of a of second- 
order in h on the grid G(x) and A 2 is a centered finite-difference approximation of a 2 = d XyX 
of second-order in h. Then J(A X ) + XaA 2 is the negative of an M-matrix and approximates 
the operator J + Xad XyX with 0{h 2 ) accuracy. 

Proof For simplicity, let us consider a uniform grid with a step size h as generalization for 
the nonuniform grid is straightforward. 



1. Approximation. Rewriting Eq.(17) in terms of the operator a and taking a series 
expansion in h, we find that 



A< 



1 + fioa + aa + 



pOi 



[l-p)9 2 



1 



0i 



+ \ha 2 ) 



+ 



9 2 + {a 
(l-p)6 2 (a 



\ha 2 ) 



-h 



h( a + \ha 2 



\ha 2 ) 



+ {a-\ha 2 )f 



h-{a + \ha 2 )f 
-- J + Xad x , x + 0(h 2 ) 



regardless of what kind of second-order approximation of a 2 in h we chose. 

2. Stability and positivity of the solution. Observe that under the conditions of Eq.(15), 
Q\I — A F is an M-matrix, and therefore its inverse is a positive matrix (see Berman and 
Plemmons (1994)). Moreover, [Oil — A F )~ 2 is also a positive matrix. The same is true for 



9 2 I + A and its inverse and inverse squared. Now, according to the previous item, we have 
freedom in choosing the discrete representation of a 2 . Let us rewrite the last three terms in 



the Eq.(17) in the form 



Ms = h< 



(1-P)02 l 

Q x hI-U F 9 2 hI + U B 2 



p9i 



pOi{U 



F\2 



B\2 



hhi - u F y 



+ 



[i-p)e 2 (u B ) 
(9 2 hl + u B ) 2 
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and the first three terms in the form 

M 4 = -I + ^(U F + U B ) + ^U F U B . 

As h — } 0, the matrix M 3 is positive. That is because (— U F )~ l is upper triangular with 
all elements on the main and upper diagonals equal to 1, and (f/ B ) _1 is lower triangular with 
all elements on the main and lower diagonals equal to 1. Therefore, the diagonal elements 
of M 3 as h — > are di = [p6\ + (1 — p)6> 2 ]/2 > 0. On the other hand, M 4 is a tridiagonal 
matrix with negative elements on the main diagonal and positive elements at the first upper 
diagonal. Since the dominant elements at the first lower diagonal are proportional to 1/h 2 
and positive, it is always possible to choose small enough h such that all elements on this 
diagonal are positive, and thus the total matrix J(A X ) is the negative of an M-matrix. Then, 



based on Proposition 4.1, the stability and positivity of the solution follow. 



Using Proposition 5.2, it is easy to construct an appropriate scheme for obtaining the 
solution of the pure jump equation in Eq.([9]). As the jump operator in this equation doesn't 
contain the term a\d x>x , the trick is to borrow it from the diffusion part. For instance, in 
the diffusion equation there is a term \<J 2 d xx . Therefore, we can split the whole operator 



5.2 



2 

J + T> by leaving half of this term (i.e., \(J 2 d x ^ x ) in V, while moving the other half to the 

jump operator. Therefore, the operator in the exponent at the second splitting step becomes 

J + \a 2 d x ^ x . A comparison of this expression with the analogous term in Proposition 

immediately gives rise to the explicit representation a = cr 2 /(4A). 

Finding a good choice of the step size h* is determined by two conditions. The first one 

requires that all elements of M 4 on the first lower diagonal must be positive. This condition 

is equivalent to 

a no 

h 2 > 2h~' 

The second condition is needed to guarantee that all diagonal elements of J(A X ) are negative. 

This gives 

2a p e 1 + (l- p )6 2 

h 2 2 

The solution of both inequalities provides the value of h* . 

5.3 CGMY Model 

Computation of jump integrals under the CGMY model (also known as the KoBoL model, 
or more generally as generalized tempered stable processes (GTSPs)) was considered in 



detail in Itkin and Carr (2012) using a similar approach. GTSPs have probability densities 
symmetric in a neighborhood of the origin and exponentially decaying in the far tails. After 
this exponential softening, the small jumps keep their initial stable-like behavior, whereas 
the large jumps become exponentially tempered. The Levy measure of GTSPs is given by 



e - v L\v\ e -VR\y\ 

KV)=*L , ,i +ttl . ly<0+Afl ■ . 1+ttH ly>0, (1* 
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where v R , v L > 0, X R , Xl > and a R , a^ < 2. The last condition is necessary to provide 



y 2 {i(dy) < oo, / fi(dy) < oo. 

1 J\y\>l 



(19) 



The next proposition follows directly from Proposition 7 of Itkin and Carr (2012) 15 
Proposition 5.3 The PIDE 



d_ 



C(x,t) 



C(x + y,T)-C(x,r)-^-C(x 1 T)(ey-l) 



l*{v)dy 



is equivalent to the PDE 




Or 



(20) 



C(x,t) = (C r + £ l )C(x,t), 

C R = X R T(-a R ) {(u R - a) aR - v** + [u%* - (u R - l) a «] a} , 

a R < 2, Re(u R - a) > 0, v R > 1, 
£_ = X L T(-a L ) {(y L + a) aL - U? + [v^ - (y L + l) aL ] a} , 

a L < 2, Re(v L + a) > 0, v L > 0, 

where T is the gamma function and Re() of an operator refers to the spectrum of the dis- 
cretization of that operator. 

In special cases, this equation changes to 



£r = X R I hg(u R ) - log (v R -a) + log ( — I a 

afl = 0,R(KR-o)>0,R(i/ fl )>l, 
C L = X L < log(z/ L ) - log [y L + a) + log ( — | a 

a L = 0, R(v L + a) > 0, R(i/ £ ) > 0, 



(21) 



and 



c R = x 



R — AR 



[u R - a) log(z/R - a) - v R log(z/ R ) + a (\og(u R - 1) - 2v R coth x (l - 2v R )) 
a R = 1, ite(z/ fl - a) > 0, z/ R > 1, 

£l = Xl 



(22) 



y u L + a) log - o(l H- z/l) log 



i>£ 



? y L 



a L = 1, Re{v L + a) > 0, z/£ > 0, 



where the logarithm of the differential operator is defined in the sense of Bakas et al. (1993). 



15 In Itkin and Carr's paper, jump integrals were denned on half-infinite positive and negative domains, 
while here they are defined on the whole infinite domain. 
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There are a few ways to proceed in this case. First, one can use an extra Strang's splitting; 



instead of directly solving Eq.(20), solve it in three sweeps. At every step, only one operator, 
either Cr or Cl enters the equation. Thus, the construction of the appropriate discrete 
operator is simplified. The second approach is based on the observation that eigenvalues of 
the sum of two M-matrices are also positive. This result follows from Wayl's inequality (see 



Bellman (I960)). Therefore, if every operator in the right-hand side of Eq.(20) is represented 



by the negative of an M-matrix, the sum of those operators is also the negative of an M- 
matrix. However, the discretization of these operators, while on the same grid, could differ, 
thus adding some flexibility to the construction of the numerical scheme. 



As shown in Itkin and Carr (2012), the computation of the matrix exponential could 
be fully eliminated by using the following approach. First, they show that for aj G Z the 
solution of the pure jump equation could be reduced to the solution of a system of linear 
equations where the matrix in the left-hand side of the system is banded. Therefore, the 
complexity of this solution is O(N). Then to compute the matrix exponential for a real 
a, first choose three closest values of aj G Z. Given the solutions at these a/, we can 
interpolate them to give the solution for a. Therefore, the complexity of this solution is that 
of interpolation. 

This approach, however, does not work well if < a < 2, since we do not have a solution 
at a = 2. To proceed in such a way would then require extrapolation instead of interpolation. 
It is well known that extrapolation is not a reliable procedure, and so in what follows we 
apply the general approach of this paper to the GTSP models. 



5.3.1 Case 1: 1< oc R < 2. 

Proposition 5.4 If 1 < a R < 2, then the discrete counterpart Lr of the operator Cr is the 
negative of an M-matrix if 



L R = \RT(-aR){(vRl-A B ) 



B\ a R 



V 



R 



+ [,%*- ( VR -1)°«]A F } 



The matrix Lr is an 0(h) approximation of the operator Cr. 

Proof In this case, r(— o,r) > 0. The last three terms in the curly braces form the negative 



of an M-matrix. Now, due to the existence conditions in Eq.(20), the step size h must be 
chosen such that Vr > 1/h, and therefore (vr — A B ) R is a nonnegative matrix. Therefore, 
if h also obeys the condition 



(y R - 1//0 



U R 



[v%« - (u R - in i/h < o, 



then Lr is the negative of an M-matrix. But this inequality is always true for h > 0. That the 
approximation is first order follows from the definition of A F and A B . Therefore, according 



to Proposition |4.1[ the above scheme is unconditionally stable and preserves the positivity 
of the solution. ■ 
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5.3.2 Case 2: a R < 0. 

Proposition 5.5 If a R < 0, then the discrete counterpart L R of the operator C R is the 
negative of an M-matrix if 

L R = X R T(-a R ) { (v R I - A B ) aR - i£* + [i/«« - (u R - l) a «] A B ] . 

The matrix L R is an 0(h) approximation of the operator L R . 



Proof The proof looks the same as that of Proposition 5A_ if one takes into account that 
V T - {vr ~ l T R < if or < 0. ■ 

5.3.3 Case 3: < a R < 1. 

The most difficult case is when < a R < 1. This is because Y{oc R ) < in this range of a R , 



so construction of the finite-difference scheme based on Proposition |4.1| experiences some 
problems. Therefore, we must consider another approach, which is closely related to the 



concept of an "eventually positive matrix"; see Noutsos and Tsatsomeros (2008). Below we 



reproduce some definitions from this paper necessary for our further analysis. 
Definition An N x N matrix A = {a^} is called 

• eventually positive, denoted by A > 0, if there exists a positive integer k such that 
A k > for all k > k ; we denote the smallest such positive integer by k = k (A) and 
refer to ko(A) as the power index of A; 

• exponentially positive if for all t > 0, e tA = Ylh=o t ~&~ > 0> 

• eventually exponentially positive if there exists to £ [0, oo) such that e tA > for all 
t > t . We denote the smallest such nonnegative number by to = ^o(^4) and refer to it 
to (A) s the exponential index of A. 



We also need the following Lemma from Noutsos and Tsatsomeros (2008): 
Lemma 5.6 Let A £ ]R 7VxiV . The following are equivalent: 

1. A is eventually exponentially positive. 

2. A + bl is eventually positive for some b > 0. 

3. A T + bl is eventually positive for some b > 0. 



Based on the above definitions and Lemma 5.6, we can prove the following proposition. 



Proposition 5.7 Suppose < a R < 1 and consider the following discrete approximation of 
the operator L R : 

L R = X R T{-a R ) { [v R I - A c ) aR - iff + [iff - (u R - l) aR ] A B } . 

Then the matrix L R (i) is an O(h) approximation of the operator C R , (ii) has all negative 
eigenvalues and (Hi) is eventually exponentially positive. 



Proof We prove each statement separately. 

Proof of (i): This follows from the fact that A c approximates the operator a to second 
order in h, while B does the same to first order. 

Proof of (%%): As r(— cxr) < 0, the matrix 

Mi = X R T(-a R ) {-v%* + [v a R * - {v R - \Y«] A B ) 

is the negative of an M-matrix, with negative elements on the main diagonal and positive 
elements on the first lower diagonal, while all other elements are zero. The term M 2 = 
XrT(— cur) {vrI — A c ) R evaluates to a dense matrix with negative elements on the main 
and first lower diagonal and positive elements on the first upper diagonal. The condition 
vr > \i(A c ), i = 1 . . . ,N is always valid since R(A») = for i = 1, . . . , N. 

Denote the elements on the first lower diagonal of M\ as <i 2 Li(M 1 ), % = 1, . . ., N — 1. 
Note that d*l 1 (M 2 ) < while dLi(Mi) > 0. Therefore, dLi(Af 2 + Mi) < d}_ x (M x ). For 
the elements on the main diagonal (we mark them with index 0) the opposite is true; i.e., 
d l (M2 + Mi) > do (Mi). Also observe that for small enough h, we have Mi oc 1/h while 
M-2 oc l/h R and < olr < 1. Hence Lr could be made strictly diagonally dominant by a 
proper choice of h. Moreover, all diagonal elements of Lr are negative. Therefore, according 



to Gershgorin's circle theorem (Golub and Van Loan (1983)), all eigenvalues of Lr have a 
negative real part. 

Proof of (Hi): Observe that Lr has the following properties for sufficiently small values 
oih: 

1. Elements on the main diagonal are negative. 

2. Elements on the first lower and upper diagonals are positive. 

3. Other elements are small in absolute values as compared with that on the above three 
diagonals and can have various signs. 



To prove the statement of (iii), use Lemma 5.6 to choose b large enough such that elements 



on the main diagonal of A + bl are all positive. Then d l (A + bl) > 0, d l _ 1 (A + bl) > and 
d\(A + bl)>0 for all i = 1, . . . , N. Also, 

|(A + 6/) 4i |< \(A + bl) u \, \fj >i + l, or Vj < i - 1, (23) 

\{A + bI)iA < |(A + 6/) M _i|,Vj >i+l,orVj<i-l. 

Taking the square of A + bl propagates large positive values on the diagonals do, d±, oLi 
to the diagonals ofo, d-2, so the elements on these diagonals become positive. That is because 



of the above inequalities in Eq.(23). Therefore, (A + bl) 1 is a positive matrix. Thus, by 
Lemma |5.6[ A is eventually exponentially positive. ■ 

Our experiments show that one can choose Ar = h sufficiently small, as in the proof of 



Proposition 5.7 to get e ArLR > for all At > Atq. Strictly speaking, the above scheme 
is not unconditionally stable because of this condition on At. However, in contrast to, 
for example, explicit finite difference schemes, this condition is not restrictive. Indeed, it 
only sets the lower limit for At, while for an explicit FD scheme the convergence condition 
determines the upper boundary of the time step. 
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5.3.4 Case 4: a R = 0. 



When a R = 0, the CGMY model reduces to the VG model (Madan et al. (1998)). The jump 



operator Jl + Jr is given in Eq.(|2T]). For the following, we rewrite the jump operator as 

vr-\* 



Cr = \ 



r 



log 1 



a 
vr 



+ a log 



VR 



When vr is not too close to 1, it makes sense to move the term A^log ( ^— ^ J a to 

the diffusion part T>r, and then approximate it as J\ = \r log ( VR ~ X ) A B (to first order 

approximation in h). This guarantees that if D is an appropriate discrete approximation of 
T> (e.g., D is the negative of an M-matrix), then D + J\ is also the negative of an M-matrix 

(because log f ^— - j < 0, and thus J± is also the negative of an M-matrix). 

Once we move the first term in Eq.(21 ) into the diffusion operator, the second (remaining) 



term in Eq.(21) could be directly exponentiated to give 



ArJR,: 



V R 



vr! 



AtAj 



A discrete approximation of this term to first order in h is 

AtAjj 

/ i / i • \ 

Q 



At l r;. 



vr 



v r I - A F 



(24) 



As {vrI — A F ) is an M-matrix, jz^f is a positive upper triangular matrix with all the 
diagonal elements satisfying Xa < 1, i = 1, . . . , N. Accordingly, Q is an upper triangular 
matrix satisfying the same property for its diagonal elements. Thus, the spectral norm of Q 
is IIQH < 1. Therefore, the proposed approximation is unconditionally stable and preserves 
positivity of the solution. 

In the case when vr is close to 1, moving J\ into the diffusion discrete operator D could 
result in needing a very small h to make the discrete diffusion operator stable. Otherwise, 



one will face a convection-dominated problem that requires special treatment; see, e.g., Duffy 



(2006). Therefore, we consider another approach below. 



Let us again take the operator J R 2 , and recall that A# log 



vr 



is an upper triangular 



u R I-AF t 

positive matrix. By construction, each nondiagonal element is less than the main diagonal 
element. This inspires the following representation. Set B 

\ogB = \og[b hl (I + K)], 



vr 



u R I-A F 



SC 



1(1 



where K is an upper triangular matrix with zero elements on the main diagonal and all other 



elements kij < 1, i < j, i,j 



N. Then 



log B = log(6i,i)7 



T ^ K 2 K 3 

2 3 



Hi 



Again for simplicity, we consider only a uniform grid here. 
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Due to the last condition on the elements of K, this series converges to a positive upper 
triangular matrix with zeros on the main diagonal. Accordingly, since < &i 7 i < 1, the 
matrix log-B is the negative of an M-matrix. This immediately implies that the discrete 



operator Jr = B+Ji is the negative of an M-matrix. Therefore, according to Proposition |4.1 
the above scheme is unconditionally stable and preserves positivity of the solution. 
We have proven the following proposition. 

Proposition 5.8 Suppose cvr = 0. Let us discretize the operator Cr as 



jr 



A, 



V R~ l \ aB , ,__ ( V R 



log -= A" + lo, 



vr J & \v R I-A F 



Then Lr is an 0(h) approximation of the operator Cr, while the whole scheme is uncondi- 
tionally stable and preserves positivity of the solution. 

5.3.5 Approximations of Cl 

Approximations to Cl can be constructed in a way similar to those corresponding to Cr. 
Below we will present a few propositions that specify our construction. Proofs of these 
propositions are omitted because they are very similar to that for Cr. 

Proposition 5.9 If 1 < oll < 2, then the discrete counterpart Ll of the operator Cl is the 
negative of an M-matrix if 

L L = X L T(-a L ) { {y L I + A F )° L - u? + \v? - (y L + l)^} A B } . 

Ll is an 0(h) approximation of the operator Cl- 

Proposition 5.10 If oil < 0, then the discrete counterpart Ll of the operator Cl is the 
negative of an M-matrix if 

L L = \ L T(-a L ) { {u L I + A F f L - v? + [u? - (v L + l) aL ] A F ) . 

Ll is an 0(h) approximation of the operator Cr. 

Proposition 5.11 Suppose < a^ < 1 and consider the following discrete approximation 
of the operator Cl: 

L l = \ L T(-a L ) { (u L I + A C )« L - u a L L + \v a L L - (v L + l) aL ] A F ] . 

Then t matrix Ll (i) is an 0(h) approximation of the operator Cl, (ii) has all negative 
eigenvalues and (Hi) is eventually exponentially positive. 

Proposition 5.12 If oll = 0, then the discrete counterpart Ll of the operator Cl is the 
negative of an M-matrix if 



A, 



^L + l\ lF , , ( VL 



log -^— \A" + lo 



vl ) *\v L I + A B 

Ll is an 0(h) approximation of the operator Cl 
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5.3.6 Second-Order Approximation in h 

The second-order approximations of Cr and Cl could be constructed in a manner similar to 
the one described in Proposition 5J2 To give an example, consider the case when a^ = 0. 

Proposition 5.13 Suppose that the following discretization scheme is in order: 



L R + (3A 2 = X R { log 






\B 



A" + log 



v R 



h- 



M A 



F\2 



[vrI-AF 2 (u R I-A F ) 2 



(25) 



where (3 = 
mates Cr + (3d. 



I^nlogl^ 



Then Lr + f3A 2 is the negative of an M-matrix and approxi- 



to 0(h 2 



Proof The proof is also similar to that of the Proposition 5.2 



1. Approximation. Rewriting Eq.(25) in terms of the operator a and expanding it into a 
series on h, we find that 



L R + (3A 2 = Aj 



log 



Vr- 1 



a + log 



Vr 



+ Pd XtX + 0(h 2 



vr ) ' \VrI -a, 

2. Positivity and stability of the solution. Since (vrI — A F ) is an M-matrix, 



VR 



positive matrix. As 



v R I-A" 



is a 



-h- 



vr{A 



F\2 



1 



2 {u R I-A F Y 
the diagonal elements of the matrix 



h- 



vr(U 



F\2 



H 



vr 



2 (vRhI-U F ) 2 ' 
v R {A F ? 



1 



arc 



v R I - A F 
\i{H) = hv R 



h- 



2 (u R I-A F f 
vr + 1/2 

(VR + I) r 



Thus, a value of h can be found (not too small, but not too large as well) such that H is 



a positive upper-triangular diagonally dominant matrix. Then log if is the negative of an 



M-matrix (see Proposition 5.8). 



Thus, Lr + f3A 2 is the negative of an M-matrix as well. 



As in Proposition 5.2, we can borrow the term (3d XjX from the diffusion operator when 
doing Strang's splitting. Again, suppose that in the diffusion equation there is a term \a 2 d XjX . 
Therefore, we can construct splitting of the whole operator Cr + T> by leaving under T> just 
a portion r\ of this term, i.e., \rja 2 d X)X , and moving the other portion, |(1 — rj)a 2 d XjX , to the 
jump operator. Therefore, the whole operator in the exponent at the second step of splitting 
becomes Cr + |(1 — rj)a 2 d x ^ x . Comparison of this expression with the analogous term in 



Proposition |5.13| gives rise to the explicit representation 

Vr 



V 



1 — h—z- log 



a z 



Vr- 1 



17 This strongly depends on how close Vr is to 1. Our experiments show that for vr 
size h — 0.95 is fine, while for vr — 1.01 this should be changed to h=l. 
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Test 


T 


K 


r 


q 


C/P 


£ 


p 


K 





A 


A*J 


^j 


1 


1 


100 


0.05 


0.0 


c 


0.3 


-0.5 


1.5 


0.1 


5 


0.3 


0.1 



Table 1: Initial parameters used in test calculations. 



6 Numerical Examples 

In the first test we used our numerical approach to compute prices of European vanilla options 
under the Bates model (a jump-diffusion model with Merton's jumps). This solution was 
compared with the semi-analytical solution obtained by using an inverse Fourier Transform 
(FFT) since the characteristic function for the Bates model is known in closed form; see, 
e.g. 



Crepey (2000). 



For the diffusion step we used the method described in detail in Hout and Foulon (2010) 



A nonuniform space grid was constructed in both x and v dimensions which contained 100 
nodes in x G [0,S max ], S max = 40max(S , , K), and 40 nodes in v G [0,v max ], v max = 5v . 
Here K is the strike, So, t>o are the initial levels of the stock price and instantaneous variance. 



For the jump step this grid was extended to S, 



up 



10 4 



Further increase of S up does not 



influence the option price much, so this boundary was chosen based on a practical argument. 
The steps of the jump grid when outside of the diffusion grid (where they both coincide with 
each other) grew according to geometric progression hi = hx g\ where h = (S max — S m i n )/N 
is an average step size for the diffusion grid, g is the growth factor, which in our experiments 
was chosen as g = 1.03. The total jump grid thus contained 237 nodes, 75 of which were the 
diffusion grid nodes. 

The initial parameters used in the test are given in Table [I] 

We computed European option prices under the Bates model in two ways. The first 
approach utilizes the fact that the characteristic function of the Bates model is known in 
closed form. Therefore, pricing of European options can be done using any FFT algorithm. 



Here we used a standard version of the Carr and Madan (1999) method with a constant 
dumping factor a = 1.25 and N = 8192 nodes. The second approach (FDE) uses an 
algorithm described in this paper, i.e., splitting and matrix exponentials, where the diffusion 



(Heston) equation was solved using the method of fractional steps described in Hout and 
Foulon| ( [20"l0"| . 

In Fig. [T] absolute and relative differences in prices obtained in our experiments are 
presented as a function of moneyness M = Sq/K. It is seen that FDE provides a very 
reasonable accuracy compared to the almost exact solution obtained with the FFT method. 

To see how much of the observed numerical error could be attributed to the Heston model 
itself, e.g., to the FD algorithm for computing a pure diffusion part, we repeated this test 
with no jumps and presented these results in Fig. [2} 

In the second test we considered a model similar to Bates, but with jumps simulated 
using the VG model. We used the parameters in Table [TJ In addition, the VG model 
parameters were chosen as: 9 = 0.1, a — 0.4, v = 3, which translates 18 to vr = 1.5098, vl = 



IN 



For explicit formulae to provide this translation, see Madan et al. (1989) 
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Figure 1: Absolute and relative differences in 
call option price as a function of moneyness 
M for the Bates model computed using an 
FFT algorithm (FFT) and the algorithm of 
this paper (FDE). 
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Figure 2: Absolute and relative differences in 
call option price as a function of moneyness 
M for the Heston model computed using an 
FFT algorithm and FDE. 



2.7598, Xr = Xl = 0.33. The grid was constructed as it was in the previous test. However 
the upper boundary of the jump grid was moved to 10 5 , and S max = 20m&x(S ,K). Again 
we computed European option prices in two ways. As the characteristic function of the VG 
model is known in closed form, the characteristic function of our model is a product of that 
for the Heston and VG models. We then used an FFT algorithm proposed by Alan Lewis, 
and as applied to the VG model discussed in detail in Itkin (2005). The second approach 
uses the FDE algorithm described in this paper. 

In Fig. [3| the absolute and relative differences in prices obtained by these two methods 
are presented as a function of the moneyness M = Sq/K. Here FDE behaves worse than 
in the case of the Bates model, because we used just the first order approximation in h. 
Still, the relative difference with the exact solution is less than 0.5%, and for M ~ 0.5 the 
difference rises to only 1.7%. 



7 Conclusion 



In this paper (which is a further extension of our paper Itkin and Carr (2012)) we proposed 
a new method to solve jump-diffusion PIDEs. This method exploits a number of ideas, 
namely: 

1. First, we transform a linear non-local integro-differential operator (jump operator) 
into a local nonlinear (fractional) differential operator. Thus, the whole jump-diffusion 
operator J + T> is represented as a sum of the linear and non-linear parts. 
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Heston + VG model 




Figure 3: Absolute and relative differences in call option price as a function of moneyness 
M under the Heston+VG model computed using Lewis's FFT algorithm and FDE. 



Second, operator splitting on financial processe^j is applied to this operator, namely 
splitting a space operator into diffusion and jumps parts. For nonlinear operators, 
this approach was elaborated on based on the definition of Lie derivative (see Koch 



and Thalhammer (2011)). The described splitting scheme provides a second-order 



approximation of J + T> in time. 

3. At the third step various finite-difference approximations of the non-linear differential 
operator J are proposed for the Merton, Kou and GTSP (a.k.a., CGMY or KoBoL) 
models. We demonstrated how to construct these approximations to (i) be uncon- 
ditionally stable, (ii) be of first- and second-order approximation in the space grid 
step size h and (iii) preserve positivity of the solution. The results are presented as 
propositions, and the corresponding proofs are given based on modern matrix analy- 
sis, including a theory of M-matrices, Metzler matrices and eventually exponentially 
positive matrices. 

Numerical examples in the paper illustrate good efficiency and accuracy of this method. 
Moreover, the performance of the method is not worse than that of FFT (see the discussion 
in the introduction). 

All these results seem to be new. Also, to the best of our knowledge, all the approaches 
to solving jump- diffusion PIDE known in the literature were either 0(h) or 0(t), while the 
proposed methods include finite-difference schemes of the second order approximation, i.e., 
0(h 2 ) + 0(t 2 ). The method is naturally applicable to both uniform and nonuniform grids, 

19 This is similar to splitting on physical processes, e.g., convection and diffusion, which is well-known in 
computational physics. 
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and is easy for programming. The method also does not require final reinterpolation (as in 
the FFT method), because the jump grid includes the PDE grid as a subset. 

Also notice that the present approach allows pricing some exotic, e.g., barrier options as 
well. In addition, it respects not just European (vanilla) payoffs but also digital, American 
and Bermudan ones. 
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